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The near-wake vortex flowfield from a NACA0012 half-wing was simulated using a fully 
unstructured Navier-Stokes flow solver in three dimensions at a chord Reynolds number of 
4.6 million and a Mach number of approximately 0.15. Several simulations were performed 
to examine the effects of boundary conditions, mesh resolution, and turbulence scheme on 
the formation of wingtip vortex and its downstream propagation. The standard Spalart- 
Allmaras turbulence model was compared with the Dacles-Mariani and Spalart-Shur 
corrections for rotation and curvature effects. The simulation results were evaluated using 
the data from an experiment performed at NASA Ames’ 32 in x 48m low speed wind tunnel. 
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Nomenclature 

half wingspan 

initial vortex pair separation 
chord 

coefficient of static pressure 
lift coefficient 
vortex circulation 
initial vortex descent velocity 
reference temperature 
reference freestream velocity 
Reynolds number 
Mach number 
vortex core radius 
streamwise direction 
vertical direction 
span wise direction 
axial velocity 


I. Introduction 

W AKE vortex studies conducted in the past can be divided roughly into two focus areas: near- wake and far- 
wake. The near- wake region extends about 10-50 wingspans behind the generating aircraft and is dominated 
by the roll-up of multiple vortices off the wing surface, flaps, and tail. At greater distances, i.e., in the far- wake 
region, these multiple vortices coalesce into a counter-rotating vortex pair and the ambient meteorological 
conditions dominate vortex transport and decay. The development of operational tools and concepts of operations 
has been driven by the far- wake research. These studies have relied on either large eddy simulations (e.g., Han et al. 
2000; Proctor 2009) or field experiments (Hinton 1995; Perry et al. 1997) and have resulted in a wealth of 
knowledge leading to a better understanding of wake decay and transport under varying conditions of atmospheric 
turbulence and stratification. The experimental and simulation data have been utilized for the development of fast- 
time wake prediction models (Sarpkaya et al. 2001; Holzapfel 2006; Proctor and Hamilton 2009). The fast-time 
wake models can be used for the systems level design of advanced air traffic management (ATM) concepts that 
safely increase the capacity of the National Airspace System (NAS). It is also envisioned that at some later stage of 
maturity, these models could be used operationally, not only within the terminal airspace but also as onboard tools to 
support concepts such as dynamic separation of aircraft. 
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Near- wake behavior is important for at least two reasons: 1) understanding and deriving initial conditions for far- 
wake models, and 2) assessing the hazard of a near- wake encounter such as due to a blunder from a parallel 
approach. The near-wake studies have focused on understanding vortex roll-up and have been limited to a small 
range of flow conditions and Reynolds numbers. Research on how to use the information from near- wake analyses 
in operational tools and concepts has also been limited. One reason for this restriction can be attributed to the high 
computational cost of running near-wake simulations. In recent years, great strides have been made in high- 
performance computing and computational fluid dynamics (CFD) techniques. The availability of powerful 
computing resources allows researchers to attempt the detailed analysis of wake turbulence in the near field of an 
aircraft. The CFD analysis can not only lead to a better understanding of the vortex roll-up process, but can also 
provide better initial conditions for the large eddy simulation (LES) codes and the fast-time wake models - thus 
bridging the gap between the near- wake research and the far- wake research efforts. For example, fast-time wake 
models require as inputs the initial vortex pair separation (bo) and the initial vortex pair descent velocity (Vo), which 
are usually based on the assumption of elliptic wing loading (Betz 1933). The large eddy simulation codes are also 
initialized with simple idealized models, based on a fully rolled-up vortex pair. The radial distributions of tangential 
velocity (v 0 ) and circulation (r) in these idealized models are represented by various relationships. The idealized 
wake vortex models use vortex core radius (r c ) as one of the inputs and various estimates ranging from 1% to 7% of 
the wingspan have been suggested in literature (Delisi et al. 2003). Even less understood are the appropriate axial 
velocity profiles for LES initialization. In the wind tunnel studies (Chow et al. 1997) as well as field experiments 
(Figure 1), significant axial velocity components have been measured and observed. The CFD simulations can 
provide information on the intensity and distribution of the axial velocity behind the aircraft. The lack of flow field 
and wake behavior data at early stages can introduce uncertainties in the design of ATM concepts such as closely 
spaced parallel runway (CSPR) operations. In blunder scenarios, e.g., during CSPR operations, an understanding of 
wakes which are not fully rolled up would be essential for determining the safety envelope of the procedure. 

Several experimental studies have been conducted in the past to characterize wake rollup from the wingtip 
(Chigier and Corsiglia 1972; Takahashi 1987; Chow et al. 1997; Gingras et al. 2001; Delisi 2012). These studies 
have provided not only insight into the wake phenomenon but also useful data for the validation of numerical codes. 
The main limitation of these studies is that they have been conducted at a small range of relatively low Reynolds 
numbers and using few test configurations. During actual flight conditions, however, the Reynolds numbers are 
much higher than can be obtained in experiments. Previous near- wake numerical studies (Dacles-Mariani et al. 
1995; Churchfield et al. 2011) have focused on the effects of mesh resolution, boundary conditions and turbulence 
schemes on the accuracy of the computed solution of the roll-up from a wingtip. In some recent studies, attempts 
have been made to simulate the wake rollup using the complete aircraft. For example, Takashi et al. (2012) used 
fuselage- wing geometry to generate initial conditions for LES and Jurkovich (2011) simulated the wake behind a 
KC-135 tanker. In the present study, the National Aeronautics and Space Administration’s (NASA) Fully 
Unstructured Navier-Stokes 3D (FUN3D) code (Anderson and Bonhaus 1994; Anderson et al. 1995; Nielsen 1998) 
is used to simulate in detail the wake in the near field of a NACA0012 half- wing. In the following sections, 
NASA’s FUN3D code and the wind tunnel test are briefly described, and the results of numerical simulations are 
presented. 



Figure 1. Wake vortex experiment conducted at NASA Langley Research Center in the early 1970s. Axial 
flows can be seen at the center of the vortex. 
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II. The Fully Unstructured Navier-Stokes 3D (FUN3D) Code 

The Fully Unstructured Navier-Stokes 3D is a versatile code that has been used to solve complex flow problems 
ranging from hypersonics to aeroacoustics (http://fun3d.larc.nasa.gov/index.html). FUN3D is based on a node- 
centered, finite-volume scheme designed to compute both compressible and incompressible flows on grids with 
different types of elements (tetrahedral, pyramid, prism and hexahedral elements, mixed-elements). The solver uses 
an artificial compressibility method for incompressible flows (Anderson et al. 1995). Several schemes are available 
for computing the convective fluxes. Computations of viscous fluxes on tetrahedral meshes are based on the Green- 
Gauss theorem and calculations on non-tetrahedral meshes are based on an edge-based approach to avoid odd-even 
decoupling. A number of turbulence models have been implemented in the code that range from single equation 
models to hybrid Reynolds -Averaged Navier-Stokes (RANS)/LES capability. Time integration is based on a 
backwards -Euler scheme with local time-stepping to accelerate convergence. 

III. Wind Tunnel Test 

Chow, Zilliac and Bradshaw (1997) conducted detailed wind tunnel studies to characterize the turbulent flow 
field in the near wake of a half- wing. The experiments were conducted in the 32 in x 48 in low speed wind tunnel at 
the Fluid Mechanics Laboratory of NASA Ames Research Center. A rectangular half- wing model with NACA 0012 
section was used in the study (Figure 2). The half- wing had a span of 3ft and a 4 ft chord (aspect ratio of 0.75). The 
wingtip was rounded (half-body of revolution) - the constant chord section of the wing section was 33.12m. The 
maximum extension of the wingtip near the quarter chord was 2.88 in. The wing trailing edge was squared-off 
(0.121m thickness). The angle of attack of the half- wing was set to 10° relative to freestream. The quarter chord 
point of the wing model was 19m downstream of the test section inlet. The freestream velocity, u ^ was 170 ft/s 
(51.82 m/s) and the Reynolds number based on the chord length was 4.6 million. The flow was made fully turbulent 
by placing a l/8m-wide line of 0.017m-diameter roughness elements on the upper and lower wing surface at a 
distance of 2 in from the leading edge. Detailed measurements were taken using static pressure taps on the wing 
surface. Turbulence, velocity, and pressure measurements were taken at various planes along the wind tunnel 
ranging from x/c = -0.591 to x/c = 0.678, where c is the wing chord length. Triple hot-wire anemometry was used 
for turbulence measurements in the wake region and a higher density of probes was placed in the vortex core region. 

IV. Results 


A. Simulation Setup 

Dacles-Mariani et al. (1993; 1995; 1996) conducted a series of simulations for this case to support the wind 
tunnel experiment. Their computational domain consisted of the entire test section of the wind tunnel. The length 
of the domain was 1.813c (87m), the width was 1.0c (48m), and the height was 0.667c (32 in). The simulation setup 
suggested by Dacles-Mariani et al. (1995) placed the inflow and outflow boundaries very close to the wing. A 
similar setup was used by Nichols (2008), Nichols et al. (2008) and Uzun et al. (2005). To avoid problems at the 
inflow and outflow boundaries, some recent studies (Kim and Rhee 2005; Churchfield and Blaisdell 2009; 
Churchfield and Blaisdell 2011; Goparaju and Zhuang 2012) have recommended larger domains. In this study, the 
domain was similar to the one suggested by Churchfield and Blaisdell (2011) in which the inflow boundary was 
placed at 2.5c from the trailing edge and the outflow boundary at 3.0c from the trailing edge (Figure 2). The meshes 
were generated using the VGRID software (Pirzadeh 1996). The viscous spacing on the wing was set to lxl0‘ 6 c. 
The tunnel walls had the same viscous spacing. A maximum stretch ratio of 1.1 was specified in the viscous region 
for the generation of volume grid. In previous studies, the viscous spacing on the wind tunnel walls was much 
coarser compared to the spacing on the wing surface. 

Earlier studies have emphasized the need for higher-order schemes for resolving the vortex core. Dacles- 
Mariani et al. (1995) recommended using at least a 5 th -order spatially accurate scheme. Nichols et al. (2008) 
compared two 5 th -order weighted essentially nonoscillatory (WENO) schemes with a 3 rd -order scheme and found 
that the 3 rd -order scheme was too diffusive and not able to accurately resolve the wingtip vortex and its propagation. 
The mesh spacing in the vortex core region in these studies was set to 0.01c in the stream wise direction and 0.005c 
(0.24m) in the span wise direction. This placed 11 to 20 points across the vortex core depending on how the mesh 
refinement was prescribed in the vortex core region. In the current study, a high mesh resolution was specified in 
the vortex core region by placing cylindrical refinement sources along the vortex core line. Two meshes (Figure 3) 
were generated with resolutions of 0.2m and 0.1 in around the vortex core. 


3 

American Institute of Aeronautics and Astronautics 


The inflow Mach number based on the freestream conditions was 0.1509 and the reference temperature, T ref was 
set to 20°C. The Reynolds number was 4.6 million. Total temperature and total pressure based on the wind tunnel 
data were specified at the inflow boundary and static pressure was specified at the outflow boundary. Although the 
Mach number is low, several studies in the past have used compressible solvers without pre-conditioning (Nichols et 
al. 2008; Churchfield and Blaisdell 2009) and obtained results comparable to incompressible solutions (Dacles- 
Mariani et al. 1995). FUN3D can be used for both compressible and incompressible flows and both options were 
tested. The results of the compressible solver are presented in this paper. The solution convergence was based on 
reduction in the residuals by several orders of magnitude and forces on the wing to become steady. The axial 
velocity was also monitored at a reference point inside the tunnel (at this point, u/u ^ = 1 in the wind tunnel data). 
The reference point is at v = -1.135c, y = 0.342c and z = 1.05 b. Convergence was reached within 3500-5000 
iterations for all cases. 

Dacles-Mariani et al. (1995; 1996) found that the standard turbulence models were not able to accurately predict 
the vortex pressure and velocity as well as the vortex core radius. The schemes predicted diffused vortices with 
increasing core radii. They suggested modifications be made to the production term of the turbulence model to 
decrease excessive diffusion in the vortex core. The Dacles-Mariani modification has since been used in several 
studies with favorable results. Kim and Rhee (2005) evaluated several turbulence schemes and found that the 
Reynolds Stress Transport Model (RSTM) described in Kim (2001) and the Spalart-Allmaras (SA) model with 
Dacles-Mariani correction gave the best results. Uzun et al. (2005) used LES for this case but at Re = 0.5 million, 
and therefore their computational results did not agree well with the experimental data. In this study, the SA model 
is examined with both Dacles-Mariani et al. (1995; 1996) and Shur et al. (2000) corrections for rotation and 
curvature effects. 




Figure 2. The computational domain for the wing -wind-tunnel configuration is shown. The inflow and 
outflow boundaries were extended in the simulations. 


Figure 3. Mesh with coarse resolution (0.2m) around the vortex core near the wingtip is shown in the left 
panel and the mesh with fine resolution (0.1m) around the vortex core region is shown in the right panel. 

B. Comparison with Wind Tunnel Data 
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The wind tunnel coordinate systems are used in the plots presented in this section. In the wind-tunnel 
experiment, the z-axis is spanwise, the y-axis is vertical, and the v-axis is in the streamwise direction. 

Surface Static Pressure Coefficient 

Comparison of the computed C p with the wind tunnel data at different locations over the half-wing is shown in 
Figure 4, and the contours of the computed C p over both the upper and lower wing surfaces are shown in Figure 5. 
The solutions are for the fine mesh using the Dacles-Mariani correction in the SA model. The simulated results are 
in good agreement with the experimental data. However, near the wingtip, there are discrepancies due to the 
inadequacy of the model to fully resolve the vortex formation. Figure 6 shows the comparison of computed C p 
using the standard SA model with both the Dacles-Mariani and the Spalart-Shur corrections. Near the wing root, all 
three models gave identical results but near the wingtip the C p prediction was slightly improved by using the 
Spalart-Shur correction. 








Figure 4. Comparison of the computed static pressure coefficient (C p ) with wind tunnel data at different 
locations over the half-wing. Solution on the fine mesh using SA model with Dacles-Mariani correction. 



I 



Figure 5. Computed coefficient of static pressure (C p ). The left panel shows the contours on the suction side 
(upper surface) of the wing and the right panel shows the contours on the pressure side of the wing. Solution 
on the fine mesh using SA model with Dacles-Mariani correction. 
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The computed lift coefficient (C L ) on the fine mesh using Dacles-Mariani correction was 0.548, compared to 
0.549 reported by Dacles-Mariani (1996) and 0.84 reported in Goparaju and Zhuang (2012). The experimental C L 
was 0.51, but Chow et al. (1997) have expressed reservations on the accuracy of the measured value given the wind 
tunnel test setup. In the wind tunnel experiment, the flow was tripped and made fully turbulent by placing 
roughness elements on the upper and lower wing surfaces near the leading edge. Dacles-Mariani (1996) included 
transition in their computations but reported no difference in results and therefore computed only the fully turbulent 
flow. Rumsey and Lee-Rausch (2012) accounted for transition which resulted in an improved C L prediction. In this 
study, fully turbulent flow was assumed and the effect of including transition is something that can be looked into in 
a future study. 




Figure 6. Comparison of the computed static pressure coefficient (C p ) with wind tunnel data near the wing 
root (left panel) and near the wingtip (right panel). In the legend, T stands for solution on the fine mesh. SA 
is the standard Spalart-Allmaras model, SADM is the SA model with Dacles-Mariani correction and SARC is 
the SA model with Spalart-Shur correction. 


Axial Velocity 

Figure 7 shows cutplanes of computed axial velocity normalized by u ^ at different locations. The computed axial 
velocity is compared with measurements along the vortex core line in Figure 9. The sensitivity to mesh resolution 
and the turbulence scheme is clearly observed in the computed solutions. The coarse mesh was not able to resolve 
the initial formation of the vortex at the wingtip and, in the case of standard SA model, the vortex dissipates rapidly 
- even on the fine mesh. The predicted values of the peak axial velocity at the trailing edge for the fine mesh were 
1.61^00, 1. 76woo, and for the SA, SA-Dacles-Mariani, and the SA-Spalart-Shur models, respectively. The 

maximum measured axial velocity was 1.7 lu^ near the trailing edge of the wing. 
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Figure 7. Cutplanes of the computed axial velocity normalized by u at different locations. Solution on the 
fine mesh (SA model with Spalart-Shur correction). 
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Vortex Static Pressure Coefficient 

Figure 8 shows cutplanes of the computed coefficient of static pressure at different locations. The computed C p 
is compared with measurements along the vortex core line in Figure 9. Use of a coarse mesh and the standard SA 
model (even on the fine mesh) resulted in large errors in the predicted values of C p . The SA model with both 
Dacles-Mariani and Spalart-Shur corrections on the fine mesh are in good agreement with measurements until the 
trailing edge but, after the vortex lifts off from the wing, the Spalart-Shur correction gives better values for C p . The 
predicted values of the minimum C p at the trailing edge for the fine mesh were -2.95, -3.3, and -3.4 for the SA, SA- 
Dacles-Mariani, and the SA-Spalart-Shur models respectively. 
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Figure 8. Cutplanes of static pressure coefficient (C p ) at different locations. Solution on the fine mesh (SA 
model with Spalart-Shur correction). 



Figure 9. Comparison of the computed static pressure coefficient (C p ) and the normalized axial velocity 
(m/«oo) with wind tunnel data along the vortex core centerline. In the legend, T stands for the fine mesh and 
6 c’ stands for the coarse mesh. SA is the standard Spalart-Allmaras model, SADM is the SA model with 
Dacles-Mariani correction, and SARC is the SA model with Spalart-Shur correction. 
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Vortex Core Location 

The predicted vortex core locations are compared with measurements in Figure 10. The predicted positions by 
different models were similar and in reasonable agreement with the experimental data. The SA model with Spalart- 
Shur correction on the fine mesh gave a slightly better prediction of the vortex core location. 




Figure 10. Comparison of the computed vortex core locations with the wind tunnel data. In the legend, 6 P 
stands for the fine mesh and 6 c’ stands for the coarse mesh. SA is the standard Spalart-Allmaras model, 
SADM is the SA model with Dacles-Mariani correction, and SARC is the SA model with Spalart-Shur 
correction. 


V. Summary 

Numerical simulations of the near-wake flowfield were performed using a fully-unstructured Navier-Stokes flow 
solver, and the results were compared against wind tunnel measurements. The solution was sensitive to boundary 
conditions and the inflow and outflow boundaries had to be placed farther out to avoid numerical problems. A mesh 
refinement study showed that a 0.0021c grid size resolution was needed in the vortex core region to adequately 
resolve the formation of the vortex at the wingtip and its downstream propagation. This resolution is more than a 
factor of 2 higher than previously suggested values (Dacles-Mariani et al. 1996). Dacles-Mariani et al. (1996) used 
a fifth-order scheme whereas the current study was performed using a second -order accurate flow solver. 
Simulations were run with the standard SA model and SA with the Dacles-Mariani and Spalart-Shur corrections. 
All turbulence models were able to predict the vortex location reasonably well and there was less dependence on 
mesh resolution. However, there were substantial differences in the prediction of axial velocity and C p in the vortex 
core. The SA model with Spalart-Shur correction gave the best predictions in the wake region and also improved 
the C p near the wingtip. The modified term in the Dacles-Mariani correction had almost no effect in the boundary 
layer but improved the solution in the wake region. Future extension of this work will focus on the simulation of 
wingtip vortex using either the Common Research Model (CRM) geometry or the Drag Prediction Workshop 
(DPW) geometry in order to provide initial conditions for large eddy simulations. 
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